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^fotract-The  proper  design  of  RF  pulses  in  magnetic  resonance 
imaging  has  a  direct  impact  on  the  quality  of  acquired  images. 
Several  techniques  have  been  proposed  to  obtain  the  RF  pulse 
envelope  given  the  desired  slice  profile.  Unfortunately,  these 
techniques  do  not  take  into  account  the  limitations  of  practical 
implementation  such  as  limited  amplitude  resolution.  Moreover, 
implementing  constraints  for  special  RF  pulses  on  most  techniques 
is  not  possible.  In  this  work,  we  propose  an  approach  for  designing 
optimal  RF  under  theoretically  any  constraints.  The  new 
technique  poses  the  RF  pulse  design  problem  as  a  combinatorial 
optimization  problem  and  uses  efficient  techniques  from  this  area 
to  solve  this  problem.  In  particular,  an  objective  function  is 
proposed  as  the  norm  of  the  difference  between  the  desired  profile 
and  the  one  obtained  from  solving  the  Bloch  equations  for  the 
current  RF  pulse  design  values.  Two  global  optimization 
techniques  were  implemented  using  genetic  algorithms  and 
simulated  annealing.  The  results  show  a  significant  improvement 
over  conventional  design  techniques  and  suggest  the  practicality  of 
using  of  the  new  technique  for  clinical  use. 

Keywords  -  MRI,  RF  pulse  design,  genetic  algorithms,  simulated 
annealing. 

I.  Introduction 

Magnetic  resonance  imaging  (MRI)  is  one  of  the  most 
promising  imaging  modalities.  It  offers  true  volumetric 
acquisition,  ability  to  visualize  and  quantify  flow,  and 
spectroscopic  imaging  to  image  both  anatomy  and  function. 
This  technique  relies  on  collecting  the  signal  from  an  excited 
slice  or  volume  within  the  human  body.  This  excitation  is 
achieved  using  special  RF  pulses  that  are  designed  to  provide 
the  required  localization  within  the  imaged  volume.  The  proper 
design  of  RF  pulses  is  important  to  avoid  artifacts  such  as 
cross-talk  in  the  acquired  images. 

The  design  of  RF  pulses  is  a  rather  difficult  problem 
mathematically  [1].  The  basic  goal  of  this  design  is  to  enable 
the  desired  slice  profile  to  be  achieved  using  the  available  RF 
pulse  generation  hardware.  The  practical  implementation  of  RF 
pulse  systems  consists  of  a  computer  that  stores  the  array  of 
digital  values  representing  the  RF  pulse  envelope  amplitudes. 
These  values  are  converted  using  a  digital-to-analog  converter 
(DAC)  into  actual  voltage  levels.  This  DAC  has  a  limited 
resolution  in  both  amplitude  and  time.  As  a  result,  the  generated 
envelope  voltages  appear  like  a  piecewise  constant  curve  with  a 
fixed  time  step  and  limited  stepwise  amplitudes.  These  levels 
modulate  the  amplitude  of  the  output  of  an  RF  generator  before 
applying  this  output  to  the  RF  coils.  The  RF  coils  may  be  either 
linear  (i.e.,  allowing  only  the  real  component  to  be  applied)  or 
quadrature  (i.e.,  allowing  both  the  real  and  imaginary 
components  to  be  used).  The  actual  frequency  of  the  RF 
generator  and  the  applied  slice  selection  magnetic  field  gradient 
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determine  the  position  of  excited  slice.  On  the  other  hand,  the 
amplitudes  of  the  RF  pulse  envelope  points  determine  the  shape 
of  the  excitation  profile  as  well  as  its  flip  angle. 

For  low  flip  angles,  the  problem  of  determining  the  RF 
pulse  profiles  can  be  approximated  by  the  Fourier  transform  of 
the  desired  slice  profile.  An  example  of  that  is  the  use  of  Sine- 
shaped  RF  pulse  envelope  to  achieve  a  gate-like  slice  profile. 
As  the  flip  angle  gets  higher,  this  approximation  becomes 
inaccurate  and  more  advanced  pulse  design  techniques  must  be 
used  [1,2].  Among  the  many  design  methodologies  that  have 
been  proposed  in  the  last  decade,  the  Shinnar-Le  Roux  (SLR) 
method  is  probably  the  most  widely  used  [2].  This  method 
works  by  transforming  the  problem  into  one  of  designing  a 
finite  impulse  response  (FIR)  filter.  The  solution  of  this 
problem  is  obtained  as  an  FIR  filter  coefficients  and 
subsequently  transformed  back  into  the  desired  RF  pulse 
envelope. 

Given  the  extensive  literature  on  FIR  filter  design,  the 
strength  of  the  SLR  technique  is  that  it  taps  into  some  of  the 
most  powerful  FIR  filter  design  techniques  to  solve  the  original 
problem  of  RF  pulse  design  [2].  Nevertheless,  it  still  has  some 
limitations  that  did  not  enable  its  performance  to  be  optimal  in 
practice.  In  particular,  this  technique  does  not  take  into 
consideration  the  limited  amplitude  resolution  in  designing  the 
RF  filter.  As  it  is  well  known  in  FIR  filter  design  literature,  the 
limited  precision  (i.e.,  quantization  or  limited  word  length 
effects)  in  implementing  the  digital  filter  may  substantially 
deteriorate  its  performance  [3].  Even  if  the  FIR  filter  design 
technique  is  modified  to  takes  care  of  this  problem  and  provide 
optimal  filter  coefficients  for  a  given  precision,  there  is  no 
guarantee  that  the  backward  transformation  to  RF  pulse 
coefficients  would  preserve  this  property  for  RF  pulse 
coefficients.  In  other  words,  the  characteristics  of  the  SLR 
transformation  do  not  allow  such  constraints  to  be  imposed.  In 
fact,  it  is  generally  difficult  to  impose  any  type  of  constraints  on 
the  solution  (like  for  example  adiabatic  constraints).  As  a  result, 
the  obtained  design  may  in  fact  be  suboptimal  in  many  cases 
that  are  common  in  practical  use.  An  example  where  difficulty 
to  obtain  accurate  slice  profiles  is  reported  is  the  use  of 
unconventional  spatial  encoding  techniques  such  as  wavelet 
encoding  and  pseudo-Fourier  imaging  [4].  The  implementation 
of  such  techniques  had  to  compromise  between  the  need  to  use 
low  flip  angles  to  obtain  accurate  slice  profiles  for  correct 
encoding  and  the  need  for  high  flip  angles  for  better  signal-to- 
noise  ratio.  Therefore,  a  new  RF  pulse  design  technique  that 
can  incorporate  practical  constraints  thus  offering  a  true  optimal 
performance  under  the  practical  implementation  constraints 
would  be  rather  helpful  to  solve  these  problems. 
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In  this  work,  we  formulate  the  problem  of  RF  pulse  design 
as  a  combinatorial  optimization  problem  with  an  arbitrary 
number  of  constraints.  This  formulation  takes  into  account  the 
limited  precision  of  RF  pulse  generation  and  provides  the 
optimal  results  at  any  given  precision.  Unlike  SLR  technique, 
the  objective  function  used  is  based  on  the  computed  slice 
profile,  which  offers  a  feedback  loop  to  improve  the  results. 
Two  optimization  techniques  are  applied  to  solve  this  problem, 
namely,  genetic  algorithms  and  simulated  annealing.  We 
provide  the  detailed  implementation  details  for  each  and  present 
their  results  compared  to  those  of  the  SLR  technique. 

II.  Theory 

Given  the  definition  of  the  RF  pulse,  it  is  possible  to 
compute  the  expected  slice  profile  using  the  solution  to  the 
Bloch  equations.  This  solution  relies  on  using  the  analytical 
form  for  the  slice  profile  from  a  single  rectangular  pulse  of 
arbitrary  magnitude  given  in  [1].  The  Bloch  equations  relates 
the  derivative  of  the  magnetization  M  to  the  current 
magnetization  value  in  the  form, 

'  0  ~G-_  By  ' 

M=y-  Gz  0  -Bx  M  =  A-M.  (1) 
-Bv  Bx  0 

Here,  G:  is  the  slice  selection  gradient,  Bx  and  By  are  the  two 
quadrature  components  of  the  RF  pulse,  y  is  the  gyromagnetic 
ratio.  For  a  rectangular  pulse,  the  solution  can  be  simply 
computed  as, 

M (z)  =  M 0(z)  ■  Exp(-A  •  t)  ,  (2) 

where  t  is  the  duration  of  the  rectangular  RF  pulse  and  the 
matrix  exponent  can  be  computed  analytically  for  this  problem 
as  in  [1]. 

Keeping  in  mind  the  practical  implementation  of  RF  pulses 
in  the  form  of  piecewise  constant  envelope  pulses  (i.e.,  a 
sequence  of  rectangular  pulses  of  arbitrary  amplitudes),  the 
output  magnetization  from  one  piece  serves  as  the  initial 
condition  for  the  next.  Hence,  given  any  design  for  the  RF 
pulse,  the  slice  profile  can  be  computed  this  method.  Given  that 
the  amplitudes  of  the  RF  pulses  must  be  represented  within  a 
certain  number  of  bits,  the  problem  now  becomes  the  one  of 
finding  the  optimal  combination  of  amplitudes  that  would  give 
a  slice  profile  closest  to  the  desired.  This  problem  description 
shows  that  this  problem  is  indeed  a  combinatorial  optimization 
problem.  Using  the  rich  literature  of  this  area,  the  solution  can 
be  obtained  efficiently  and  accurately.  In  this  work,  we  explore 
two  of  the  most  prominent  techniques  in  this  area,  namely, 
genetic  algorithms  and  simulated  annealing. 

III.  Methods 

A.  Optimization  Using  Genetic  Algorithms  (GA) 

The  GA  method  employs  mechanisms  analogous  to  those 
involved  in  natural  selection  to  conduct  a  search  through  a 
given  parameter  space  for  the  global  optimum  of  some 
objective  function  [5].  The  main  features  of  this  approach  are: 


•  A  point  in  the  search  space  is  encoded  as  a  chromosome. 

•  A  population  of  N  chromosomes/search  points  is 
maintained. 

•  New  points  are  generated  by  probabilistically  combining 
existing  solutions. 

•  Optimal  solutions  are  evolved  by  iteratively  producing  new 
generations  of  chromosomes  using  a  selective  breeding 
strategy  based  on  relative  values  of  the  objective  function  for 
the  different  members  of  the  population. 

A  solution  is  encoded  as  a  string  of  genes  to  form  a 
chromosome  representing  an  individual.  In  many  applications 
the  gene  values  are  [0,1]  and  the  chromosomes  are  simply  bit 
strings.  An  objective  function,  f,  is  supplied  which  can  decode 
the  chromosome  and  assign  a  fitness  value  to  the  individual  the 
chromosome  represents. 

Given  a  population  of  chromosomes  the  genetic  operators 
crossover  and  mutation  can  be  applied  in  order  to  propagate 
variation  within  the  population.  Crossover  takes  two  parent 
chromosomes,  cuts  them  at  some  random  gene/bit  position  and 
recombines  the  opposing  sections  to  create  two  children.  For 
example,  for  one-point  crossover  crossing  the  chromosomes 
010-11010  and  100-00101  at  randomly  selected  position  3-4 
gives  010-00101  and  100-1 1010.  Another  example  is  for  two- 
point  crossover,  where  crossing  the  chromosomes  01-01 1-010 
and  10-000-101  at  randomly  selected  two  positions  2-3  and  5-6 
gives  01-000-010  and  10-011-101.  Mutation  is  a  background 
operator,  which  selects  a  gene  at  random  on  a  given  individual 
and  mutates  the  value  for  that  gene  (for  bit  strings  the  bit  is 
complemented). 

The  search  for  an  optimal  solution  starts  with  a  randomly 
generated  population  of  chromosomes  and  an  iterative 
procedure  is  used  to  conduct  the  search.  For  each  iteration,  a 
process  of  selection  from  the  current  generation  of 
chromosomes  is  followed  by  applying  the  genetic  operators  and 
re-evaluation  of  the  resulting  chromosomes.  Selection  allocates 
a  number  of  trials  to  each  individual  according  to  its  relative 
fitness  value.  The  fitter  an  individual  the  more  trials  it  will  be 
allocated  and  vice  versa.  Average  individuals  are  allocated  only 
one  trial.  Trials  are  conducted  by  applying  the  genetic 
operators  (in  particular  crossover)  to  selected  individuals,  thus 
producing  a  new  generation  of  chromosomes.  The  algorithm 
progresses  by  allocating,  at  each  iteration,  ever  more  trials  to 
the  high  performance  areas  of  the  search  space  under  the 
assumption  that  these  areas  are  associated  with  short 
subsections  of  chromosomes  which  can  be  recombined  using 
the  random  cut-and-mix  of  crossover  to  generate  even  better 
solutions  [5]. 

The  use  of  GA  is  robust  in  that  they  are  not  affected  by 
spurious  local  optima  in  the  solution  space.  Nevertheless,  The 
parameters  that  control  the  GA  can  significantly  affect  its 
performance,  and  there  is  no  guidance  in  theory  as  to  how 
properly  select  them.  The  most  important  parameters  are  the 
population  size,  the  crossover  rate,  the  mutation  probability  and 
the  fitness  function.  The  following  are  the  proposed  algorithm 
details. 


1)  Population  of  the  chromosomes'.  Population  represents  the 
size  of  the  solutions  that  we  are  working  with.  The  bigger 
is  the  population  size,  the  faster  the  convergence,  but  the 
more  the  computational  burden  at  the  initial  iterations.  It 
can  be  fixed  or  dynamically  changed  allover  the  run  of  the 
algorithm.  In  our  problem,  a  fixed  population  size  of  1000 
individual  has  been  used. 

2)  Population  Initialization :  The  chromosomes  are  initialized 
randomly.  A  chromosome  that  represents  the  SLR  solution 
is  added  to  the  initial  population. 

3)  Chromosome  structure'.  Binary  chromosome  is  used.  The 
RF  pulse  is  encoded  into  the  chromosome  as  follows.  The 
real  RF  pulse  values  are  converted  into  discrete  ones 
according  to  the  resolution  of  the  D/A  of  the  MRI  machine 
(12  or  16  bit  for  example).  Every  bit  represents  a  gene.  The 
most  significant  bits  of  all  values  are  placed  adjacent  to 
each  other,  then  the  second  most  bits  and  so  on  until 
placing  the  least  significant  bits  together  at  the  end.  Flence, 
the  chromosome  size  equals  the  number  of  the  RF  pulse 
envelope  values  times  the  bit  resolution  (usually  12  or  16). 

4)  Fitness  criterion'.  The  chromosome  is  decode  to  obtain  the 
corresponding  RF  pulse  envelope  amplitudes  and  its  slice 
profile  is  computed  by  solving  the  Bloch  equations.  Then, 
an  error  measure  is  calculated  for  the  difference  between 
the  response  of  this  RF  pulse  and  the  desired  response. 
This  measure  is  usually  taken  as  either  the  1 -norm  or  2- 
norm  of  the  difference  vector.  The  results  of  this  paper 
were  obtained  using  the  1-norm. 

5)  Selection  scheme'.  A  biased  roulette-wheel  selection  is 
used.  In  this  method,  the  fitness  of  all  the  chromosomes  in 
the  population  is  evaluated.  Divide  each  of  these  fitness 
values  by  the  total  summation  of  these  fitness  values  and 
then  multiply  it  by  100  to  get  the  percentage  that  this 
chromosome  can  be  evolved  in  the  next  generation. 

6)  Crossover.  One  point  crossover  is  used  the  probability  of 
crossover  is  taken  as  90%. 

7)  Mutation'.  The  probability  of  mutation  is  taken  as  small  as 

1%. 

B.  Optimization  Using  Simulated  Annealing  (SA) 

The  concepts  of  simulated  annealing  are  based  on  a  strong 
analogy  between  the  physical  annealing  process  of  solids  and 
the  problem  of  solving  large  combinatorial  optimization 
problems.  Annealing  is  a  thermal  process  for  obtaining  low 
energy  states  of  a  solid  in  a  heat  bath.  It  contains  two  steps. 
First,  the  temperature  of  the  heat  bath  is  increased  to  a 
maximum  value  at  which  the  solid  melts.  Second,  it  is  decrease 
carefully  until  the  particles  arrange  themselves  in  the  ground 
state  of  the  solid.  In  the  ground  state  the  particles  are  arranged 
in  a  highly  structured  lattice  and  the  energy  of  the  system  is 
minimal.  The  ground  state  of  the  solid  is  obtained  only  if  the 
maximum  temperature  is  sufficiently  high  and  the  cooling  is 
done  sufficiently  slow.  Otherwise  the  solid  will  be  frozen  into  a 
meta-stable  state  rather  than  ground  state. 

The  Metropolis  algorithm  is  a  numerical  technique  that 
simulates  the  annealing  process  [6].  In  this  algorithm,  given  a 


current  state  i  of  the  solid  with  energy  £j,  then  a  subsequent 
state  j  is  generated  by  applying  a  perturbation  mechanism 
which  transforms  the  current  state  into  a  next  state  by  a  small 
distortion  (for  instance  by  displacement  of  a  particle).  The 
energy  of  the  next  state  is  Er  If  the  energy  difference  is  less 
than  or  equal  to  0,  the  state  j  is  accepted  as  the  current  state.  If 
the  energy  difference  is  greater  than  0,  the  state  j  is  accepted 
with  a  certain  probability  which  is  given  by  Exp((£j— EJ)l 
( Kb.T )),  where  T  denotes  the  temperature  of  the  heat  bath  and 
Kb  is  a  physical  constant  known  as  the  Boltzmann  constant.  The 
acceptance  rule  described  above  is  known  as  the  Metropolis 
criterion  and  the  algorithm  that  goes  with  it  is  known  as  the 
Metropolis  algorithm. 

We  can  apply  the  Metropolis  algorithm  to  generate  a 
sequence  of  solutions  of  a  combinatorial  optimization  problem. 
For  this  purpose,  we  assume  an  analogy  between  a  physical 
many-particle  systems  and  a  combinatorial  optimization 
problem.  This  is  based  on  the  equivalencies  that  solutions  in  the 
combinatorial  optimization  problem  are  equivalent  to  states  of  a 
physical  system,  and  that  the  cost  of  a  solution  is  equivalent  to 
the  energy  of  a  state. 

The  control  parameter  here  is  equivalent  to  the  temperature. 
That  is,  if  the  energy  difference  is  greater  than  0,  the  state  j  is 
accepted  with  a  certain  probability  which  is  given  by  the 
Metropolis  criterion  where  KB.  T  is  replaced  by  a  single  control 
parameter  c.  Initially  the  control  parameter  contains  a  large 
value.  This  means  that  most  changes  even  those  with  large 
deterioration  to  the  objective  function  will  be  accepted.  The 
value  of  c  is  gradually  decreased  with  each  step  until  it 
approaches  a  very  small  value  near  0.  Then  no  deterioration  to 
the  objective  function  will  be  possible.  This  feature  allows  the 
simulated  annealing  algorithm  to  escape  from  local  minima  [6], 

In  order  to  apply  SA  to  our  RF  pulse  problem,  we  start  with 
the  RF  pulse  designed  by  SLR  Algorithm.  A  new  solution  is 
generated  by  adding  a  vector  of  random  numbers  to  the  vector 
containing  the  RF  pulse  values.  The  random  vector  should  be 
with  small  values  to  generate  a  new  solution  with  energy  (cost 
value)  with  small  distortion  of  the  old  solution.  The  cost 
function  used  is  the  same  as  the  evaluation  function  used  in  the 
GA.  The  control  parameter  c  is  decreased  linearly  with  each 
step.  The  constraints  on  the  solution  are  applied  to  the  new 
values  in  every  step  to  make  sure  that  the  solution  is  acceptable 
(e.g.,  the  constraints  of  limited  amplitude  resolution).  The  new 
values  are  accepted  or  rejected  according  to  the  outcomes  from 
the  objective  function  and  the  Metropolis  criterion.  The  process 
in  repeated  until  the  changes  in  the  solution  become  very  small. 

VI.  Results  and  Discussion 

The  proposed  techniques  were  applied  to  design  two  RF 
pulses  with  rectangular  spatial  profiles  at  ji/2  and  71  flip  angles. 
The  outcome  of  the  SLR  technique  is  compared  to  the 
outcomes  from  optimized  RF  pulses  using  the  GA  and  SA 
methods.  The  results  are  shown  in  Figs.l  and  2.  The  1-norm  of 
the  error  between  the  desired  and  obtained  is  shown  in  Fig.  3. 
As  can  be  observed,  both  techniques  show  a  significant 
improvement  over  the  design  computed  using  the  SLR 


technique.  Also,  the  design  using  the  SA  method  appears  better 
than  that  from  the  GA  technique.  Even  though  both  techniques 
should  theoretically  reach  the  global  optimum,  we  have  to 
realize  that  both  techniques  are  iterative.  This  means  that  the 
user  decides  when  the  iteration  stops  both  in  time  and  accuracy. 
Hence,  the  above  results  from  comparing  GA  and  SA 
techniques  indicate  that  the  SA  technique  reaches  the  solution 
faster  than  GA  in  the  number  of  experiments  we  performed. 

One  of  the  most  important  applications  for  the  proposed 
technique  is  in  the  area  of  unconventional  spatial  encoding 
where  complex  RF  pulse  are  to  be  generate  at  high  accuracy. 
The  proposed  technique  is  expected  to  enable  better 
reconstruction  accuracy,  less  image  artifacts,  and  higher  signal- 
to-noise  ratio.  The  study  of  this  area  requires  the  investigation 
of  several  of  these  techniques  and  the  assessment  of  the  results 
and  is  left  to  future  work. 

The  addition  of  constraints  on  the  solution  using  the 
proposed  approaches  is  straightforward.  In  each  of  the 
algorithms  for  GA  or  SA,  the  obtained  new  solution  after 
perturbation  is  passed  through  a  ‘filter’  that  checks  the  validity 
of  the  new  solution  from  the  point  of  view  of  the  required 
constraints.  In  very  much  the  same  way  this  was  used  for 
practical  amplitude  resolution,  other  conditions  can  also  be 
included  for  specialized  RF  pulse  such  as  adiabatic  pulses.  The 
investigation  of  this  application  should  be  addressed  in  future 
work. 


V.  Conclusions 

A  new  optimized  RF  pulse  design  approach  is  proposed. 
The  new  technique  relies  on  posing  the  problem  as  a 
combinatorial  optimization  problem  and  uses  genetic 
algorithms  and  simulated  annealing  to  compute  the  solution 
under  any  type  of  constraints.  The  results  demonstrate  the 
success  of  the  new  approach  and  suggest  its  potential  for 
practical  use  in  clinical  magnetic  resonance  imaging. 
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Figure  1 :  Comparison  between  slice  profiles  for  flip  angle  of  90  degrees. 
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Figure  2:  Comparison  between  RF  pulses  for  flip  angle  of  180  degrees. 
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Figure  3:  Comparison  between  error  in  slice  profiles  versus  flip  angle. 


